NewGridFloatFromNetCDF Subroutine

private subroutine NewGridFloatFromNetCDF(layer, fileName, variable, stdName, time)

Uses

create a new raster_real reading data from NetCDF file The variable to read can be defined by its current name or the standard_name. The dimensions x and y of the variable is calculated searching from the dimensions of the couple of variables with 'standard_name' equal to 'projection_x_coordinate' and 'projection_y_coordinate' for projected reference systems or 'longitude' and 'latitude' for geographic reference systems or 'grid_longitude' and 'grid_latitude' for rotated pole systems If a comprehensible reference systems is not found a geodetic reference system is supposed. Once the variable is retrieved, offset and scale factor are applied and a check on minimum and maximum valid value is performed.

Arguments

Type IntentOptional Attributes Name
type(grid_real), intent(out) :: layer

gridreal to return
character(len=*), intent(in) :: fileName

NetCDF file to be read
character(len=*), intent(in), optional :: variable

variable to read
character(len=*), intent(in), optional :: stdName

standard name of the variable to read
type(DateTime), intent(in), optional :: time

time of the grid to read

Variables

Type Visibility Attributes Name Initial
character(len=80), public :: attribute
integer(kind=short), public :: i

loop index
integer(kind=short), public :: ios

error return code
integer(kind=short), public :: j

loop index
integer(kind=short), public :: latlon
integer(kind=short), public :: nDimsVar

number of dimensions of a variable
integer(kind=short), public :: nVars

number of variables
integer(kind=short), public :: ncId

NetCdf Id for the file
integer(kind=short), public :: ncStatus

error code returned by NetCDF routines
real(kind=float), public :: offset
real(kind=float), public :: scale_factor
character(len=1), public :: shp
integer, public, ALLOCATABLE :: slice(:)
real(kind=float), public, ALLOCATABLE :: tempGrid(:,:)
integer(kind=short), public :: varId

variable Id
character(len=100), public :: variableName

Source Code

SUBROUTINE NewGridFloatFromNetCDF &
!
(layer, fileName, variable, stdName, time)

USE StringManipulation, ONLY: &
!Imported routines:
StringCompact


IMPLICIT NONE

! Arguments with intent(in):
CHARACTER (LEN = *), INTENT(in) :: fileName  !!NetCDF file to be read
CHARACTER (LEN = *), OPTIONAL, INTENT(in) :: variable  !!variable  to read
CHARACTER (LEN = *), OPTIONAL, INTENT(in) :: stdName  !!standard name of 
                                                      !!the variable  to read
TYPE (DateTime), OPTIONAL, INTENT(in) :: time  !!time of the grid to read

! Arguments with intent(out):
TYPE (grid_real), INTENT (out)            :: layer  !!gridreal to return


! Local scalars:
INTEGER (KIND = short)          :: ios !!error return code
INTEGER (KIND = short)          :: ncStatus !!error code returned by NetCDF routines
INTEGER (KIND = short)          :: ncId  !!NetCdf Id for the file
INTEGER (KIND = short)          :: varId !!variable Id
INTEGER (KIND = short)          :: nVars !!number of variables
CHARACTER (LEN = 80)            :: attribute
INTEGER (KIND = short)          :: nDimsVar !!number of dimensions of a variable
INTEGER (KIND = short)          :: i, j   !!loop index
CHARACTER (LEN = 100)           :: variableName
CHARACTER (LEN = 1)             :: shp
REAL (KIND = float)             :: scale_factor
REAL (KIND = float)             :: offset
INTEGER (KIND = short)          :: latlon

! Local arrays:
INTEGER, ALLOCATABLE            :: slice (:)
REAL (KIND = float), ALLOCATABLE :: tempGrid (:,:) 

!------------end of declaration------------------------------------------------

!------------------------------------------------------------------------------
![1.0] open NetCDF dataset with read-only access
!------------------------------------------------------------------------------
ncStatus = nf90_open (fileName, NF90_NOWRITE, ncId)
IF (ncStatus /= nf90_noerr) THEN
  CALL Catch ('error', 'GridLib',        &
  TRIM (nf90_strerror (ncStatus) )//': ',  &
  code = ncIOError, argument = fileName )
ENDIF

!------------------------------------------------------------------------------
![2.0] get x and y size and allocate array
!------------------------------------------------------------------------------

CALL GetXYSizesFromFile (ncId, layer % jdim, layer % idim, shape = shp, latlon = latlon)

!allocate grid
ALLOCATE ( layer % mat (layer % idim, layer % jdim), STAT = ios )
IF (ios /= 0) THEN
  CALL Catch ('error', 'GridLib',  &
  'memory allocation ',  &
  code = memAllocError,argument = variable )
ENDIF


!allocate temporary grid
IF (latlon == 1) THEN !coordinate order lon-lat
  ALLOCATE ( tempGrid (layer % jdim, layer % idim), STAT = ios )
ELSE !coordinate order lat-lon
  ALLOCATE ( tempGrid (layer % idim, layer % jdim), STAT = ios )
END IF

IF (ios /= 0) THEN
  CALL Catch ('error', 'GridLib',  &
  'memory allocation ',  &
  code = memAllocError,argument = variable )
ENDIF

!------------------------------------------------------------------------------
![3.0] get values
!------------------------------------------------------------------------------

IF ( PRESENT (variable) ) THEN
  ncStatus = nf90_inq_varid (ncId, variable, varId) 
  CALL ncErrorHandler (ncStatus)
ELSE IF (PRESENT (stdName) ) THEN !search variable corresponding to standard name

  !inquire dataset to retrieve number of dimensions, variables 
  !and global attributes
  ncStatus = nf90_inquire(ncId,  nVariables = nVars )
                          
  CALL ncErrorHandler (ncStatus)
  
  DO i = 1, nVars
    attribute = ''
    ncStatus = nf90_get_att (ncId, varid = i, name = 'standard_name', &
                             values = attribute)
    IF (ncStatus == nf90_noerr) THEN !'standard_name' is found
      IF ( StringCompact(attribute(1:LEN_TRIM(attribute))) == stdName(1:LEN_TRIM(stdName)) )THEN
         varId = i
      END IF
    END IF
  END DO
END IF

!verify number of dimensions of variable to read
ncstatus = nf90_inquire_variable(ncId, varId, ndims = nDimsVar) 
CALL ncErrorHandler (ncStatus)

!If number of dimensions = 3 read info about time
IF (nDimsVar == 3) THEN
  ALLOCATE ( slice(3) )
  !Read time information
  CALL ParseTime (ncId, layer % reference_time, layer % time_unit)


  !Define current time slice (if present)
  IF (PRESENT (time) ) THEN
    slice = (/ 1, 1, 1 /)
     slice(3) = TimeIndex (ncId, layer % reference_time, &
                          layer % time_unit, time)
     

    layer % current_time = time
    layer % time_index = slice(3)
  ELSE !read the first grid
    slice = (/ 1, 1, 1 /)
    layer % time_index = 1
    !set current time
    layer % current_time = TimeByIndex (ncId, layer % reference_time, &
                                layer % time_unit, layer % time_index) 
  ENDIF
  !define time of the next grid
  layer % next_time = NextTime (ncId, layer % reference_time, &
                                layer % time_unit, layer % time_index)                 
ELSE
   ALLOCATE ( slice(2) )
   slice = (/ 1, 1/)
   layer % current_time = timeDefault
   layer % next_time = timeDefault
END IF 

ncStatus = nf90_get_var (ncId, varId, tempGrid , start = slice)
CALL ncErrorHandler (ncStatus)

!transpose temporary matrix to grid specification
CALL SwapGridRealForward (tempGrid, layer % mat, latlon)

!deallocate temporary grid
DEALLOCATE (tempGrid)

!------------------------------------------------------------------------------
![4.0] get attributes
!------------------------------------------------------------------------------
ncStatus = nf90_get_att (ncId, varId, name = 'standard_name', &
                             values = layer % standard_name)
                            
ncStatus = nf90_get_att (ncId, varId, name = 'long_name', &
                             values = layer % long_name)                              
                                                       
ncStatus = nf90_get_att (ncId, varId, name = 'units', &
                             values = layer % units) 
                            
ncStatus = nf90_get_att (ncId, varId, name = 'varying_mode', &
                             values = layer % varying_mode)
!If attribute is not defined set to default                             
IF (ncStatus /= nf90_noerr) layer % varying_mode = 'sequence'                            
                               
ncStatus = nf90_get_att (ncId, varId, name = '_FillValue', &
                             values = layer % nodata)

!if _FillValue is not defined search for missing_value   
IF (ncStatus /= nf90_noerr)  THEN                         
  ncStatus = nf90_get_att (ncId, varId, name = 'missing_value', &
                             values = layer % nodata)
END IF   

 !if missing_value is not defined set to default
IF (ncStatus /= nf90_noerr)  THEN                         
  layer % nodata = MISSING_DEF_REAL
END IF   
   

                         
ncStatus = nf90_get_att (ncId, varId, name = 'valid_min', &
                             values = layer % valid_min) 
!If attribute is not defined set to default                             
IF (ncStatus /= nf90_noerr) layer % valid_min = layer % nodata 
                                                         
ncStatus = nf90_get_att (ncId, varId, name = 'valid_max', &
                             values = layer % valid_max) 
!If attribute is not defined set to default                             
IF (ncStatus /= nf90_noerr) layer % valid_max = layer % nodata
                                                        
ncStatus = nf90_get_att (ncId, varId, name = 'scale_factor', &
                             values = scale_factor)
!If attribute is not defined set to default                             
IF (ncStatus /= nf90_noerr) scale_factor = 1.0
                             
ncStatus = nf90_get_att (ncId, varId, name = 'add_offset', &
                             values = offset)
!If attribute is not defined set to default                             
IF (ncStatus /= nf90_noerr)  offset = 0.  

ncStatus = nf90_get_att (ncId, varId, name = 'esri_pe_string', &
                             values = layer % esri_pe_string)
!If attribute is not defined set to default                             
IF (ncStatus /= nf90_noerr) layer % esri_pe_string = ''                         

!------------------------------------------------------------------------------
![5.0] check values and apply scale factor and offset
!------------------------------------------------------------------------------
!retrieve name of variable
ncstatus = nf90_inquire_variable(ncId, varId, name = variableName)   
layer % var_name = variableName
                          
DO i = 1, layer % idim
  DO j = 1, layer % jdim 
    IF ( layer % mat (i,j) /= layer % nodata ) THEN
    
      !apply scale factor
      layer % mat (i,j) = layer % mat (i,j) * scale_factor
      
      !add offset
      layer % mat (i,j) = layer % mat (i,j) + offset
      
      !check lower bound
      IF ( layer % valid_min /= layer % nodata .AND. &
           layer % mat (i,j) < layer % valid_min ) THEN
        layer % mat (i,j) = layer % valid_min
        CALL Catch ('info', 'GridLib', 'corrected value exceeding &
         lower bound in variable: ', argument = variableName )
      END IF  
    
      !check upper bound
      IF ( layer % valid_max /= layer % nodata .AND. &
           layer % mat (i,j) > layer % valid_max ) THEN
        layer % mat (i,j) = layer % valid_max
        CALL Catch ('info', 'GridLib', 'corrected value exceeding &
         upper bound in variable: ', argument = variableName )
      END IF  
    END IF
  END DO
END DO 

!------------------------------------------------------------------------------
![6.0] set file name
!------------------------------------------------------------------------------
layer % file_name = fileName  

!------------------------------------------------------------------------------
![7.0] read georeferencing informations from netCDF file
!      Regularize grid if necessary
!------------------------------------------------------------------------------

IF (shp == 'm') THEN !2 dimensional (matrix) coordinate for not regular grid
  IF (.NOT. ALLOCATED (layer % Iraster)) THEN
    ALLOCATE (layer % Iraster(layer%idim,layer%jdim))
    ALLOCATE (layer % Jraster(layer%idim,layer%jdim))
    CALL GetGeoreferenceFromNCdataSet (ncId, varId, layer % cellsize, &
                 layer % xllcorner, layer % yllcorner, &
                 layer % grid_mapping, layer % Iraster, layer % Jraster, .TRUE.)
  END IF
  
  
  !create a temporary copy of layer % mat
  IF (ALLOCATED (tempGrid) ) THEN
    DEALLOCATE (tempGrid)
    ALLOCATE (tempGrid(layer%idim,layer%jdim))
  ELSE
    ALLOCATE (tempGrid(layer%idim,layer%jdim))
  END IF
  tempGrid = layer % mat
  layer % mat = layer % nodata
  !fill in the regular grid
  DO i = 1, layer%idim
    DO j = 1, layer%jdim
      IF (layer % Iraster (i,j) /= -9999 .AND.  layer % Jraster (i,j) /= -9999) THEN
        layer % mat (layer%Iraster(i,j),layer%Jraster(i,j)) = tempGrid (i,j)
      ELSE
        layer % mat(i,j) = layer % nodata
      END IF
    END DO
  END DO
  DEALLOCATE (tempGrid)
ELSE !regular grid stores coordinate in 1dimensional array (vector)
  CALL GetGeoreferenceFromNCdataSet (ncId, varId, layer % cellsize, &
                 layer % xllcorner, layer % yllcorner, &
                 layer % grid_mapping, layer % Iraster, layer % Iraster, .FALSE.)                                
END IF


!------------------------------------------------------------------------------
![8.0] close NetCDF dataset
!------------------------------------------------------------------------------
ncStatus = nf90_close (ncid)
CALL ncErrorHandler (ncStatus)

END SUBROUTINE NewGridFloatFromNetCDF